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Abstract. In this article we describe applications of Discrete Differential 
Forms in computational GR. In particular we consider the initial value prob- 
lem in vacuum space-times that are spherically symmetric. The motivation 
to investigate this method is mainly its manifest coordinate independence. 
Three numerical schemes are introduced, the results of which are compared 
with the corresponding analytic solutions. The error of two schemes converges 
quadratically to zero. For one scheme the errors depend strongly on the initial 
data. 



1. Introduction 

Most methods that are presently used in numerical GR are in some sense referred 
to a coordinate system. This can be a major problem, because not only is it 
impossible in general to cover a global space-time with a single coordinate chart. 
But also it is generally impossible to know beforehand the effects that certain gauge 
conditions specified during the course of a simulation will imply. 

In view of this the question occurs, as to whether it is possible to develop a 
numerical method that is manifestly coordinate invariant. One such method is 
Regge calculus [I] which, unfortunately, so far has not played a role in computa- 
tional GR (see however [2]). In other approaches to treat the problem of coordinate 
dependencies multiple coordinate systems are used to cover the space-time [3J HI [5] . 

However, a manifestly invariant numerical method must be based on invariant 
quantities describing the geometry of space-time. The prime examples for invariant 
quantities on a manifold are the scalar fields, but in the usual description even they 
are coordinate dependent, because the description of the points of the manifold 
themselves depends on the choice of a coordinate system. Therefore, in order to 
avoid coordinates we must not even use coordinates for the localisation of points of 
the space-time manifold. This implies that we cannot use the usual definition of a 
manifold as a collection of coordinate charts with transition functions which is the 
basis of almost all analytical and numerical treatments of the Einstein equations. 

In the usual procedure for discretisation the manifold structure is untouched 
while the equations are discretised, i.e., evaluated only for a finite number of points 
of the manifold. When the use of coordinates is to be avoided one has to start 
the discretisation at an even lower level namely on that of the manifold itself. 
Hence, on a quite basic level the space-time can be considered as a collection of 
abstract objects called points. The structures on the space-time are then described 
as certain relations between these points. In our case the relevant structure consists 
of primarily topological and geometric relationships. For the present purpose we 
find it more reasonable to consider the topological relationships as given in advance 
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so that the aim of computational GR is to find the geometric relations between the 
points based on an appropriate formulation of the Einstein equations. 

To achieve this we continue the work presented in [BJ. That is we approximate 
a manifold and its differentiable structure by a cellular paving [7], i.e. a collection 
of finitely many cells. The cells are the images of a certain number of standard 
shapes like (hyper-)cubes or n-simplices. In the case where all standard shapes 
are simplices we talk about a triangulation and the cellular paving is a simplicial 
complex. The cellular paving is supposed to have the same topological properties 
as the envisaged space-time. 

To illustrate the idea we consider the example of a standard 3-simplex which 
can be viewed as the interior of a tetrahedron. It is a 3-dimensional manifold. Its 
boundary is composed of four 2-simplices (faces), six 1-simplices (edges) and four 0- 
simpliccs (nodes) . With such p-simplices we can associate several quantities which 
can be interpreted in a physical way. Examples are the charge inside a volume, 
a flux through a face, the work done along an edge or the value of a potential at 
a given point. In all these cases we associate numbers with a simplex and these 
numbers are usually obtained by integration, i.e., by adding up contributions from 
'infinitesimal' pieces making up the finite simplex. So, in each case we obtain a 
map from p-simplices to numbers. 

Differential p-forms can be viewed as 'the objects which are integrated over p- 
dimcnsional submanifolds' so they provide maps from p-dimensional submanifolds 
to the reals. Thus, the maps presented above correspond to differential forms, but 
restricted to p-simplices. These objects are known as discrete differential forms. 
They have received some attention since Bossavit [8] had pointed out that they 
correspond to the lowest order mixed finite element spaces defined by Nedelec [5] 
(see also [10]). Finite elements of mixed type have been used successfully in nu- 
merical applications to electrodynamics, see [7l[TTJ[T2]. In numerical GR the finite 
element method is used e.g. in |13j . 

Our task is now to relate geometric properties such as lengths, angles, holonomies 
and curvature using differential forms to the triangulation and the various parts of 
the simplices respectively. Since 0-forms are functions they describe properties at 
single points. In order to formulate relations between points such as the distance 
between two points or the holonomy around a loop we need p-forms with p > 0. 

In order to use this approach one needs to have a formulation of geometries and, 
in particular, of GR which uses differential forms. A formulation of geometries based 
on differential forms has been provided by E. Cartan |14j . The further step towards 
a formulation of GR using differential forms has been carried out by several authors. 
We mention here the work of Sparling [15j who has set up an exterior differential 
system of equations which is closed if and only if the vacuum Einstein equations 
hold. In [6j we have shown in detail how to set up the discrete formalism based on 
this exterior system using the ideas explained above. 

In summary, the variables of our proposed discrete formulation will be the in- 
tegrals of differential forms over submanifolds. In order to get a finite number of 
variables we use a finite number of these submanifolds based on a triangulation 
of the computational domain and discretise a description of GR that uses (finitely 
many) differential forms. 

The formulation of GR that we use is based on the Cartan formalism of moving 
frames and Sparling's exterior system for vacuum GR. In this article we describe 
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a simplification of the general formalism which occurs in spherical symmetry. The 
plan of the paper is as follows. In sect. [2] we describe the equations which result 
from a symmetry reduction. In section[3]we present three possibilities to implement 
these equations in a fully discrete evolution scheme. In section [4] we discuss how 
the method can be tested and in section [5] we present the results of those tests. 
Some final remarks can be found in section [6l 

2. The spherically symmetric equations 

We start with the formulation of GR using exterior forms [TS]. The basic vari- 
ables in this formalism are the four 1-forms of a pseudo-orthonormal tetrad 6 l , 
i = 0, ...,3 [16]. Together with the coefficients rjik — diag(l, — 1, — 1, — 1) they 
define the metric as 

(1) g = e° ® 6>° - e 1 ® e 1 - e 2 ® e 2 - e 3 ® e 3 = r/ lk e l ® e k . 

For the description of the connection in this formalism sixteen 1-forms u) l k , i,k = 
0, . . . , 3 are used. The connection should be compatible with the metric and tor- 
sion free, which translates into the antisymmetry requirement and the first Cartan 
equation, respectiveljQ: 

(2) V ik ta k j + r) jk u> k i = 0, d6 i + w' fc 0* = 0. 

Furthermore the metric should fulfil Einstein's field equations, which is equivalent 
to 

(3) OLi = Si + 87rTi*S fc . 
Here, Tn~ is the usual energy momentum tensor and 

6 

are the so called hypersurface 3-forms. The forms Lj and Si are the Nester-Witten 
2- form and the Sparling 3- form, defined by (see [B"l I15j): 

(4) U = \e ijkl ^6\ S 4 = \e m (^V m £T - ^ . 

In vacuum, when T k = 0, these equations determine the geometry of space-time. 
If there is matter, additional matter equations are needed. However, we will be 
concerned only with the vacuum case so that we will have to solve the equations 

(5a) dO l + uj l k 9 k = 0, 

(5b) dL 4 = Si. 

Although the geometry is fixed, there is still the freedom of choosing a gauge, i.e. 
there are Lorentz transformations of the tetrad that do not change the metric 

(6) g = 7 U e l ®8 k = rj lk A l J j ® = (r llk A t J A k l ) 8 3 ® 6 l . 

That means by using the tetrad for the description of geometries we introduced 
unphysical (gauge) degrees of freedom. The same problem occurs when coordinate 
systems are used. However, we believe that the tetrad, beeing a geometric object, 
has a more intuitive meaning than mere coordinates. Therefore, it might be easier 
to choose a useful tetrad than an appropriate coordinate system. 



Here and in what follows it is understood, that the product of differential forms is the anti- 
symmetrised tensor product, i.e. the exterior product. 
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In this work we will concentrate on general relativistic systems with spherical 
symmetry. Thus, we will assume that the rotation group 5*0(3) acts isometrically 
on the space-time and that the orbits of this action are 2-dimensional space-like 
submanifolds. These are necessarily spheres whose area we write as AirR 2 . In 
appendices \X\ and [Bl it is shown, how to 'factor out' the symmetry action i.e., the 
angular dependence and how to derive an exterior system on the 2-dimensional 
'orbit space' M.\ spanned by the radial and the time directions 

This can be done by a decomposition of the 4-dimensional space-time manifold 
into the 2-dimensional spheres and the 2-dimensional orbit space, followed by some 
simplifications and results in the following system 

(7a) = d/o - - ± (fl - 3/ 2 - 4 ) 8° + (/o/i) 0\ 



2 V 1 JU R 2 , 

(7b) o = dfx - f u - \ [il - 3fi + + (foh) e°, 

(7c) 0^d(v^(i? 2 /o 2 -^ 2 A 2 + l)), 

(7d) = jo (dO 1 + uG°) - A (d0° + ujO 1 ) . 

Here (O ^ 1 ) is a dyad in the 2-dimensional orbit space Mi which carries a Loren- 
tzian metric. The S0(1, 1) connection on this space is given by the 1-form oj. It 
is a consequence of the equations above that this connection is torsion free. The 
geometric properties of the orbits are described by the functions /o, fi and R 2 (e.g. 
4irR 2 is the area of the orbit) . For details see appendix [XJ 
By introducing the 1-forms 

(8) R R 

a := f a 9» + hO 1 , /3 := f^ + foO 1 

the equations (|7|) can be rewritten as follows 



(9a) 


d~8°- 


^ljO 1 +aO° = 0, 


(9b) 


dO 1 


f ojff + O.0 1 = 0, 


(9c) 




da = 0, 


(9d) 


d/3 4 


2a.(3+~e (] ~e 1 = 0, 


(9e) 


du> 


-CK/3-0V = 0, 



together with the algebraic relations 

(10) a0° + (39 l = 0, at) 1 + (30° = 0. 

The two equations (flU)) are needed to ensure that a and (3 are constructed out of 
only two functions. They can be interpreted as 

(11) *a=/3 

where * is the 2-dimensional Hodge operator [l7j. In the calculations we will use 
discrete versions of both (jTDJ) and (fTTj) . 



2 Only for _R > the orbits are 2-dimensional, which implies that the symmetry action can not 
be 'factored out' when R = 0. 
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From the definition of the Hodge operator it is clear, that *9 1 = 8°. This means 
that in ((9]) we have actually only three 1-forms a, 6 1 and uj and their Hodge duals. 
From now on the Hodge duals will be called dual forms, whereas the forms a, 
and ijj themselves get the generic name direct forms. We also have for every form 
and its Hodge dual an equation, where its exterior derivative is involved except for 
the dual of oj. It turns out, that the exterior derivative of this form is pure gauge 
because it can be set to any value by an appropriate choice of gauge in Aii- This 
property can be used to fix the gauge. 

If the geometry is regular at the origin (i? = 0) then we may also include it into 
the computational domain. Thus, the origin becomes a boundary. We choose the 
gauge such that the dyad at the origin can be defined as the limit of i.e. 
such that this limit exists. 

Then it is easy to verify that if we choose at R — an arbitrary vector V = 
V°e + V 1 e 1 with finite components (V , Vi), then di? 2 (V) = and di?(V) = V R < 
co. It follows from ([56]) that in the limit R -> the expression R(f 8° + /i^XV) 
must be finite. It is clear that i?0°(V) = and RO^V) = at R = 0. Thus, 
the limits of Rfo and Rfi must be finite there, and hence fo and /i diverge. That 
means we must use the variables fo := Rfo and f\ := Rfi instead of fo and f\. 

Now we insert this observation into the equation (|7cj) , or rather into its solution 

(12) i?(/ 2 - ft + 1) = const. 

Since fo and f\ are finite, it follows that the constant vanishes. It turns out that in 
general this constant is twice the mass of the black hole, and thus the limit R — > 
only exists for the flat geometry. In this case one may easily rewrite (|7|) and ([9|)- (fTTj) 
such that the variables are fo, fi, ot := Roc, (3 := Rf3, 9° , 1 and w. 

We are now in a position to apply the discretisation procedure as explained 
in [7]. However, the situation in GR is significantly different from electrodynamics 
so that there is no straightforward implementation. For one thing the distinction 
between direct and dual forms is used in electrodynamics in an elegant way by 
employing a dual mesh. It is not so clear whether one can make use of this also 
in GR, because the description of dual forms on the dual mesh strongly relies on 
the Euclidean character of the metric. Consequently in electrodynamics one only 
discretises space using the method of discrete differential forms while we aim at 
fully discrete schemes. 

Furthermore, electrodynamics is a linear theory while the nonlincarity of GR is 
apparent from the appearance of the wedge product which has to be implemented 
on the discrete level. So we cannot simply adapt the implementation of electro- 
dynamics, but we will have to look at the forms differently. The aim is to split a 
possibly large system of nonlinear equations into many small systems of equations 
in order to get one independent system for every simplex. How this is done in 
practice, is described in [6] and briefly in the following section. 

3. Implementation of the discrete equations 

In section [2] and appendix TK\ we derived systems of exterior equations on M\. 
Now we present methods for discretising them and develop numerical schemes. 
Since there is no unique natural way, we use different discretisation schemes to 
explore various possibilities. The first scheme will be based on the system ((7|), the 
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second one is obtained through a discretisation of the equations ([9|). (fT0|) and to get 
the third scheme we discretise the system (|9]) . (fTTj) . 

3.1. Application of Whitney forms. As indicated in the introduction, we ap- 
proximate the manifold by taking into account finitely many of its subsets. More- 
over we want to use discrete differential forms as an approximation for the contin- 
uous differential forms. 

It is known that the opposite, i.e. the extension of a discrete differential form to 
a continuous differential form can be done with the help of the so-called Whitney 
forms [T8j . These are a special class of forms which can be used to construct con- 
tinuous forms from discrete forms. However, they can exist only on special domains 
such as simplices, n-dimensional cubes and shapes, that can be constructed from a 
cube by collapsing some of its edges |19) such as pyramids or prisms. The numerical 
variables are then the integrals of the forms over the corresponding figures. 

We have chosen the first possibility, i.e. we are searching for an approximation 
of space-time by taking subsets into account, that can be continuously mapped to 
simplices. In the 2-dimensional case these simplices are nodes, edges and faces. 
The reason for this choice is, that for other shapes it is not possible to get an ex- 
terior product from (anti-)symmetry requirements alone. This leads to anaesthetic 
ambiguities, and since the exterior product of Whitney forms is in general not a 
Whitney form, symmetry assumptions are the best way to introduce the discrete 
exterior product. 

Using simplices one gets, up to a normalisation, an exterior product essentially 
from the following requirements 

(1) The discrete exterior product fulfils the usual commutation rule for forms, 
i.e. for a p-form a p and a g-form (3 q we have 

a p f\/3 q = (—l) pv /3 q A ol v . 

(2) When the orientation of the simplex is changed, the sign of the correspond- 
ing value of the discrete exterior product changes. The same is true for 
every Discrete Differential Form. 

These requirements lead almost immediately to the following formula for the dis- 
crete exterior product between 1-forms a and (3 

a.(3[n 0l n ll n 2 ] = ^ (a[n , m]/3[n , n 2 ] - a[no,n 2 ]/3[n ,ni] 
( 13 ) +a[m,n 2 ]/3[ni,no] - a[ni, n ]/3[m,n 2 ] 

+ a[n 2 ,7io]/3[n 2 ,ni] - a[n 2l n 1 }f3[n2,n a }) 1 

where the expression -y [no, . . . , n p ] is the numerical variable corresponding to the 
integral of the p-form 7 over the simplex with nodes {no, . . . , n p } and orientation 
given by the ordered tuple of vectors {n\ — no, . . . , n p — no). It turns out that this def- 
inition and its analogues for higher degree forms yield an algebraic structure which 
is not associative in contrast to the continuous case. How this non-associativity in- 
fluences the method is not clear. It is clear however, that the terms which become 
ambiguous due to the non-associativity are of higher order so they converge to zero 
faster in the continuum limit. 
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For the discretisation of the exterior derivative, we remember Stokes theorem 
and get for a 1-form a 



Figure 1. The boundary of the 2-dimensional oriented face 
[na,m,ni] is composed of three oriented edges [rio,rii], [710,712], 
[711,712] and three nodes no, ni, 712. 

3.2. Properties of the simplicial mesh. Now we come to the numerical schemes. 
Common to all three schemes is the way of generating a simplicial approximation 
of a subset of M\. We start from appropriate initial data. That means from 
somewhere we have a 1-dimensional simplicial complex Ci [20 , that approximates 
a space-like curve in M.\ (see section [4] for details). At each node of Ci two linearly 
independent light-like directions 1q and li exist. It is clear, that a light-like curve 
with tangent-vector lo at a node no will have an intersection with a light-like curve 
with tangent vector li at another node 7ii0 To create a face of the mesh, that 
contains an edge of Ci, we require the two missing edges to be approximations of 
these light-like curves. Their intersection becomes the new node 712- For obvious 
reasons, these faces are called upwards directed. 

This construction seems to be the simplest invariant method to define the po- 
sition of 71 2 in (1 + l)-dimensional manifolds. It is a geometric construction and 
can be generalised to higher dimensions. The choice of the nodes at a later time 
and their connections to the nodes at the initial time is essentially arbitrary and 
only restricted by topological considerations. It is only when the 6 l are known on 
all the edges that the geometry of the mesh is determined. We will see later that 
a part of these values can be specified freely while the rest is determined from the 
equations. 

As in all numerical simulations degeneracies may occur. For instance two adja- 
cent nodes in the same level may have a time-like distance. However, this must be 
seen as a sign that the mesh is too coarse and should be refined. 

The other type of faces is called downwards directed, and is created by joining 
the intersections of the light-like curves in adjacent upwards directed faces. When 
Ci has n edges, we now have n upwards directed and (n — 1) downwards directed 
faces. The collection of these faces will be called the first time-step. 

Obviously, this procedure can be continued by taking the collection of the non 
light-like edges of the downwards directed faces as a new initial complex C[, until 
the intersection of the light-like curves from the boundary of Ci is reached. That 

^At least when no and n\ are sufficiently close to each other. 



(14) 



da[n ,ni,n 2 ] 



a[ni,77 2 ] - a[n ,n 2 ] + a[n ,m]. 
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initial simplicial complex d 



Figure 2. The triangulation of the domain of dependence of d. 

means that we calculate the domain of dependence of Cj. In principle we could 
have implemented boundary conditions, but we wanted to concentrate on the time 
evolution scheme and periodic boundary conditions are not possible in spherical 
symmetry. Figure [2] shows the triangulation. 

Having a simplicial mesh, the exterior product and derivative, we can now take 
care of the discrete equations. 

In what follows we will present three numerical schemes. Common to all schemes 
is that for each triangle a system of equation has to be solved. These are coupled 
non-linear algebraic equations. Their analysis is somewhat complicated and their 
status is not yet clear. They might not have a unique solution. However, at least one 
solution can be found by Newton's iteration method. We used the GNU Scientific 
Library, especially the implementation of a root finding algorithm called modified 
Powell method by the developers [3TJ [55] . 

3.3. Scheme I. For the first scheme the system ([7]) is used. The variables are the 
discrete 1-forms 8°, 1 and u> as well as the discrete 0-forms /o, /i and Rr 2 . For 
the upwards directed faces the numbers 

{/oM, /oMJiMiAM' -R^M,^ 2 ^!], 0°[no,ni], ^[noim], w[no,m]} 

are given initial data, and 

{/oM, h[n 2 ], R~ 2 [n 2 ], 

6°[n 0l n 2 },9 [n 1} n 2 ], 6 1 [n Q , n 2 ], 1 [ni, n 2 ], w[n , n 2 ], w[n x , n 2 }} 

are the unknowns. 

The equation (|7c|) is the statement that the function 

(15) F := (ft - ft + RT*) / (RT*) t/a 

is constant on M.±, and it is implemented in this way. Thus, we calculate the 
constant C := F[n ] from the (known) values /ofaoL /i[«o] and R~ 2 [no] and then 
require 

(16) F[n 2 ] = ((/ 2 - fl + R- 2 ) I (R- 2 f 2 ) [n 2 ] = C, 
since equation ([7c]) implies that F[no] = F[m] = F[n 2 ] = C. 
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Therefore, the number of equations is six (two 1-form equations for the two light- 
like edges, one 2-form equation and (|16|)). but the number of unknowns is nine. We 
eliminate two unknowns by using the definition of the position of n 2 (see section 
1X21) . We get 

(17) (0° - 9 l )[n ,n 2 ] =0, {8° + 9 l )[n u n 2 ] = 0, 

expressing the fact that the edges are null. 

What remains is the freedom of choosing a gauge. This corresponds to the choice 
of a cobasis at n 2 . The dyad at n 2 can be obtained from parallel transport of 
{0°, 9 1 } from the initial hypersurface to n 2 along an edge. Since parallel transport 
is defined by lj we choose it such that 

(18) u[n ,n 2 ]=0. 

The reason to choose this gauge condition is of course that it is the most simple 
one. Clearly, in general it is the aim to choose the gauge in some sense 'optimal', 
in order to get good results. However, until now we do not understand well, what 
'optimal' means here. Probably this may become clear once the properties of the 
equations are better understood. 

In the continuum limit (|18p corresponds to a dyad that is obtained through 
parallely transporting the (known) basis at the initial hypersurface along the light- 
like curves with tangent vector (eo + ei). This fixes the so-called strong conformal 
geometry of the null-hypersurface generated by the spherical family of light rays 
(see Penrose [23]). 

For the downwards directed faces the situation is easier. The 0-form equation 
(HHJ) is automatically fulfilled at all nodes, and the 1-form equations (|7a|) and (|7bp 
are fulfilled at the light-like edges. What remains are two 1-form equations (|7aj) . 
(|7bj) and one 2-form equation (|7dj) . The unknowns are the integrals of the l-forms 
{9°, , u)} along the new edge. 

Altogether we have six equations and six unknowns for the upwards directed 
faces, as well as three equations and three unknowns for the downwards directed 
ones. That means to obtain a numerical solution we only have to solve a system 
of six equations for each upwards directed triangle and a system of three equations 
for each downwards directed one. In the simulations the root finding algorithm 
sometimes detects a solution that is a bad approximation of the analytic solution. 
However, this can be resolved by choosing other starting values for the Newton 
iteration, and we did not investigate it further. 

3.4. Scheme II. In the second scheme we use the system ([9]) together with 

(19) a0° + (3~9 l = 0, aO 1 + (39° = 0. 

In this case, the variables are the discrete l-forms a, /3, 9° , 9 and u>. 
The given initial data for an upwards directed face are 

{a[n , m],/3[n , ni], [n , ni],0 [n , m], u>[n , ni]}, 

and the unknowns 

{a [n , n 2 ] , a[ni, n 2 ], /3[n , n 2 ] , P[ni, n 2 ] , 9° [n , n 2 ] , 9° [n\ , n 2 ] , 

1 [n Q , n 2 ], 9 [m, n 2 ], u;[n , n 2 ], n 2 ]}. 
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The discretisation of the 2-form equations ([9|), (fT9|) leads to seven relations. The 
number of unknowns is ten. With the same procedure as in scheme I, i.e. using the 
fact that the new edges are light-like and choosing a gauge with oj [710,712] = 0, we 
reduce the number of unknowns to seven. 

In this case the downwards directed faces are more difficult to treat. Since all 
equations came from the discretisation of 2-forms, we still have seven equations for 
these faces, but the number of unknowns is only five (the values of five 1-forms on 
the upper edge). To get around this difficulty the following idea is used. 

In general there is no exact solution of the discrete equations, because there are 
more equations than unknownsQ However, possibly one can choose the dyad such 
that finding an exact solution is possible. To find this optimal gauge one searches 
for the exact solution and the dyad simultaneously. Hence we are using the gauge 
freedom to change the number of unknowns. 

To clarify the details of this procedure we first want to discuss which gauge 
choices can be made. As a starting point serves the discussion of the gauge issues 
in [6]. There it is argued that in the intersection of two separate regions where 
the tetrads are chosen independently exist transition maps which mediate between 
the different gauges. Obviously these transition maps have the form of (position 
dependent) Lorentz transformations. 

We may choose an open covering of the manifold with those regions, such that 
every open set in the covering contains only a single simplex. Then the transition 
maps can be interpreted as gauge transformations from one simplex to another (cf. 
figure [3|) . 



Figure 3. Each simplex is viewed as being contained in an open 
set on which a tetrad is locally defined. Transition maps medi- 
ate between different tetrad patches and hence between different 
simplex gauges. 

In fact we will use these transition maps as the new unknowns, but we will 
parameterise them through their action on the tetrad and the connection forms. 
Let the dyad and the connection in the upwards directed triangle be {f? ,^ 1 ,^;}, 
and in the downwards directed face {9° , 1 ,<*>}. Furthermore let [711,712] be the 
common edge of the two simplices. 



4 The reason for this problem is of course that equations l|19|l describing the Hodge operator 
can be localised at the faces of the mesh. This causes difficulties, because it is not natural for this 
operator. On a 2-dimensional manifold the continuous Hodge makes it possible to identify 1-forms 
with pseudo 1-forms and visa versa. Thus it would be more natural to localise this operator at 
1-dimensional submanifolds. However, having no dual mesh here, it is not clear how to achieve 
this in general. 
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In two dimensions a Lorentz transformation is completely determined by a single 
parameter, its rapidity ip. Hence the action of the transition function has the form 

(20) 8° = coshi/>0° + sinhV>6> 1 , 1 = cosh^ 1 + sinh^ , Q = u - dtp. 

On the discrete level the rapidity ip can be seen as a 0-form, i.e. a map that 
assigns a number to each node. Yet, in the intersection of two tetrad patches 
there are only two nodes ri\ and U2, and hence the discretised transition map is 
determined by two parameters ip\ = ip[ n i] an( i = V ; [ ?J 2]- 

Assuming that these parameters are small we may perform a Taylor expansion 
of (|20|) . The result is that in the leading order £? 0//1 [7ii, 712] — 9°^ 1 [ni,n2] depends 
on the sum (ipi +'02), but u>[ni, n-^ —Lj[m, ri2\ depends on the difference (ipi —ip2)- 
Hence we can choose the parameters ipi and ip2 such that the values of 6 % and l 
at [711,712] ar e the same, while the values of u) and u> differ. 

With the two differences of the gauge parameters along the two light-like edges 
of the upwards directed face we obtain two new unknowns. Effectively, this is the 
same as forgetting about the value of us at these edges. Hence the two additional 
unknowns are the values of u) at the light- like edges. 

It is known from the continuous theory that this regauging should be possible, 
but it is not known what kind of restrictions are imposed on the gauge parameter 
by the discretisation. Thus, we assume that this procedure is allowed and check 
with numerical tests whether this is true. 

3.5. Scheme III. The third scheme is based on the system (|9|) . (fTTj) . This scheme 
is a modification of scheme II in the following sense. In section [2] we discussed 
that (fT9)) can be interpreted as a = */3. When evaluated on the light-like edges of 
upwards directed faces, this formula implies 

(21) [a-0)[n o ,n a ]=O=(a + p)[ni,n a ]. 

These are equations on the two light-like edges. So, instead of two 2-form equations 
we have two 1-form equations. 

This reinterpretation of (|19p does not change the number of equations for up- 
wards directed faces, but it does change it for downwards directed ones. The new 
1-form equations are of course already satisfied at the light-like edges, so we loose 
two equations and we need no regauging, since the number of equations and un- 
knowns are already equal. 

In order to test the influence of the gauge choice at the upwards directed faces 
we make another change. We will not use ui [no, 7J2] — 0, but 

(22) u>[no,»2] + w [ni,na] = 

instead. This is an implicit and non-local definition of the gauge, but it has the 
advantage that the edges [710,712] and [711,712] are treated on a par. 

An obvious problem of this scheme is that the discretisation ((21]) of the Hodge-* 
operator can only be applied when two edges of every face are light-like. Clearly 
this limits its applicability and the question occurs whether a similar discrete Hodge 
operator exists for space-like and time-like edges. Unfortunately we did not find 
such a discretisation, but we want to point out what we believe are the first steps 
to find a more general discrete Hodge operator with similar properties as (|2ip . 

For the nodes x, y and z we may interpret the edges e = [x,y] and e = [x,z], 
as vectors in the tangent space at x. The 2-dimensional Hodge operator is defined 
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/3A*a = g- 1 (j3,a)9°e 1 V/3. 



Applying (f2l?)) to the "vectors" e and e and identifying the result of applying a and 
-ka to these vectors with a[e], a[e], -ka[e] and -ka[e] appropriately leads to 



with e l = 9 1 [e] , e 1 = 9 l [e] . This is a possible definition of a discrete Hodge operator 
that can also be applied to non light-like edges, and when e or e is light-like it 
becomes the discrete Hodge operator (|2"Tj) . 

However, we identified the edges with vectors at x. In principle there is no reason 
that e.g. e is a vector at x and not at y. But imposing the analogue of (|2"4")) at y and 
z leads to an inconsistent system. Thus the question is how to formulate a discrete 
geometry in a consistent way. This will be the topic of future investigations. 



Having clarified how we apply discrete differential forms in General Relativity, 
we now describe how we tested the obtained code. Since the spherically symmetric 
vacuum solutions of the Einstein field equations are all contained within the Kruskal 
solution [24] for some value of the mass parameter M we know the exact solutions 
for our problem. The Kruskal metric has in the standard coordinate system the 
form 

(25) g = f(R) (dT ® dT - dX <g> dX) - R 2 (d-d ® di? + sin 2 dip ® dtp) . 
There the functions f(R) and R — R(T,X) are defined through 



(26) /( J R) = ^_e-w, R(T,X) = 2M(l + W((X 2 -T 2 )/e)) 1 



where W is the Lambert W- function [5S]. To start the time evolution we need 
initial data which we obtain from one of the analytic solutions. These data must 
be given as edge- and node values. 

4.1. The continuous forms. The best way to do this is to find a description of 
the geometry, that uses (continuous) differential forms, i.e. maps from the set of 
all submanifolds to the real numbers. However, such an (abstract) map can hardly 
be useful for concrete calculations, because without coordinates it is even difficult 
to describe the position of a point. Therefore we take coordinate representations 
of the Minkowski and Schwarzschild geometries. For the Schwarzschild geometry it 
is convenient to use Kruskal coordinates, because then the light-like curves (which 
are special for the described method) take a very simple form and hence it is easier 
to compare the results. 

To obtain the differential forms, we make a gauge choice, i.e. we choose 6° and 
9 , such that they generate the corresponding metric. In Minkowski space the 
natural choice is 



(24) 




4. Test scenarios 



9° = dt, 9 1 = dr, a 



dr 



dt 



(27) 



r 



r 



w = 0, i? 2 -r 2 , f 



h 



i 



r 



where r and t are the standard space and time coordinate, respectively. 
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In Kruskal coordinates we use, with the standard space and time coordinates X 



and T, 


as well 


as the mass parameter M, 




(28a) 


e Q = 


■■ y/f(R)dT, 




= Vf(R)dx, 


(28b) 


fo = 


Th(R) 
y/WY 


h 


Xh(R) 
y/WY 


(28c) 


a = 


■ h(R){XdX — TdT), 


P 


= h{R){XdT -TdX), 


(28d) 


UJ = 


g{R){TdX - XdT), 


R 2 


= 4M 2 (l + W ({X 2 -T 2 )/e)f 



where the functions h and g are defined through 

(29a) h {R) = ™P- e -&, 

AM 2 ( 1 1 

(29b) g(R) 



R \R 2M 



e 



2M 



4.2. Getting initial values. Next an initial hypersurface has to be chosen. For 
the test we used curves, whose space-time coordinates (y ,?/ 1 ) depend linearly on 
the curve parameter. These curves will be called "straight" : 

(*)-(*M3 

for some A £ [0, 1] and fixed y§, Vqi Vi> Vi- 

We get the initial edges and nodes by subdividing this curve into pieces of equal 
'coordinate length'. That is we start at a hypersurface of the form (|30|) with bound- 
ary at A = and A = 1 and subdivide it into n pieces. The pieces are again straight, 
and A takes values in the intervals [(i — l)/n, i/n], i = 1, . . . ,n. On the edges we 
can integrate the continuous 1-forms known from the analytic solution, and the 
results are the initial values for the corresponding discrete 1-forms. 

To get initial values for the 0-forms, we evaluate the corresponding continuous 
functions at the boundary points A = i/n, i = 0, . . . ,n of the sub- intervals. These 
values are invariant under coordinate transformations so that we do not start with 
coordinate dependent values in the beginning. 

4.3. Examination of the results. In order to compare the numerical results with 
the analytical solution we need to determine the location in space-time of the nodes 
and edges used in the algorithm. This can be a difficult task because in principle 
one needs to solve the geodesic equations to obtain the light rays used to define the 
nodes in the next time-slice. However, here this is very much simplified since in 
Kruskal coordinates as well as in Minkowski coordinates the radial light rays move 
on straight lines. 

When comparing the results we need to worry about the gauge. I.e. when we 
use different gauges for the discrete approximation and for the analytic solution, 
we cannot expect to get the same results. However, if the method is feasible, 
we can expect that gauge invariant discrete variables are good approximations 
of the continuous ones. Gauge independent values are for instance the lengths 



I = y\{0 ) 2 — {O 1 ) 2 ] of the space-like edges, the values of the 1-form a. on these 
edges and the values of the function R~ 2 at the nodes. 

Another way to evaluate a numerical method is a self-convergence test. There 
one compares solutions obtained on coarse meshes with a solution that is calculated 
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on a fine mesh. It is clear that a method that converges to the analytical solution 
when the mesh is refined is self-convergent, too. 

5. Results 

For the number of initial edges we have always chosen a power of two (n — 
2 l , i = 1,2,3, . . .), and calculated half of the maximal number of time-steps: for 
n initial edges, we calculate n/2 time-steps. The last time-step then contains n/2 
downwards directed simplices. Each of these simplices contains one non light-like 
edge and each of these edges contains two nodes. Altogether these are n/2 edges 
and (n/2 + 1) nodes, since nodes of adjacent edges coincide. 

It turns out that the computational costs of the three schemes are quite com- 
parable. In scheme I the solution in a mesh with approximately 200.000 faces is 
obtained in one minute on a 700MHz PC. In schemes II and III it takes about 40% 
and 10% longer respectively. This could be expected, because in scheme I one has 
to solve systems of six and three equations for the upwards and downwards directed 
triangles respectively. In schemes II and III the sizes of these systems are seven 
and seven, respectively seven and five. Thus, in scheme I the number of equations 
is much smaller. Clearly, since there is a system of equations on each individual 
face, the time that is needed to find the solution depends linearly on the number 
of faces, and hence quadratically on the number of initial edges. 

In the left diagrams of Figs. [4] [5] and [6] we show the maxima over the n/2 edges 
of the relative errors in the values of the 1-form a for schemes II and III as well as 
the maxima over the (n/2 + 1) nodes of the relative errors in the values of R~ 2 for 
scheme 10 In the right diagrams of Figs. 0] and [5] the maxima over the n/2 edges 
of the relative error of the invariant lengths I are plotted. 
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scheme II 



10 100 1000 
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Figure 4. Maximal relative errors in Minkowski space-time. Left: 
values of R~ 2 at the nodes (scheme I) and values of a at the space- 
like edges (schemes II, III). Right: invariant lengths of space-like 
edges. 



5.1. Minkowski space-time. For the initial hypersurface we have chosen the set 
{(t, r) = (0, 1 + A) : A £ [0,1]}, with the standard time and space coordinates 



5 Since we compare relative errors it is feasible to draw different quantities in the same diagram. 
Furthermore a and i?~ 2 both are used as parameterisations of the eliminated degrees of freedoms 
in the orbits of the SO (3)-group action. 
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t and r, respectively. We compare the results with the analytic solution at the 
hypersurface {(t,r) = (0.25, 1.25 + A/2) : A € [0, 1]}. 

We see from figure 21 that the relative error converges for all three schemes 
quadratically to zero when the typical size of simplices in the mesh is decreased. 
The error of the lengths is about 100 times bigger for scheme I, but it remains small 
even for the coarsest mesh. 

5.2. Kruskal geometry. In Kruskal geometry we test the code for one space-like 
and one time-like initial hypersurface. Since in Kruskal coordinates the horizon is a 
regular null-hypersurface we can test how the code behaves near the event horizon. 
So, we take the space-like curve to cross the horizon. 

5.2.1. Space-like initial data. In (T, X)-coordinates we choose the initial hypersur- 
face {T = A/2, X = — 1 + 2A, A e [0, 1]}, and compare the results at the hypersurface 
{T = 5/8 + A/4, X = -3/8 + A, A e [0, 1]}. 




10 100 1000 
number of initial edges 

scheme I - - e - - 
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FIGURE 5. Maximal relative errors in Kruskal space-time (space- 
like initial data). Left: values of R~ 2 at the nodes (scheme I) and 
values of a at the space- like edges (schemes II, III). Right: invariant 
lengths of space-like edges. 



Fig. shows that the schemes I and III provide small errors and the relative 
error converges quadratically to zero when the size of the simplices is decreased. 
The errors in scheme II on the other hand are very big, and even become bigger 
when the mesh is refined. Clearly this points to a problem in scheme II that we 
discuss later. 

5.2.2. Time-like initial data. In (T, X)-coordinates we choose the initial hypersur- 
face {T = X,X = 3, A e [0,1]}. The comparison of the results is done at the 
hypersurface {T = 0.25 + A/2, X = 3.25, A € [0, 1]}. 

Again, as can be seen in Fig. [5J schemes I and III provide very small errors that 
converge quadratically to zero when the simplex size is reduced. With 1% for a 
and 10% for the lengths the errors of scheme II are in this case also quite small. 
However, the relative errors do not become smaller for finer meshes. 

Since in fig. [6] one sees that the errors of the results in fine meshes are nearly 
the same as the errors in coarse meshes, we also performed a self-convergence test, 
where the number of initial edges in the finest mesh was 2048. The result is that 
scheme II converges linearly to a solution that differs from the analytical one. 



16 



RONNY RICHTER, JORG FRAUENDIENER, AND MARLENE VOGEL 



Aa/a, AR~ 2 /R 
0.01 
0.001 
le-04 
le-05 
le-06 
le-07 
le-08 
le-09 
le-10 



1 1 ^ I I I ll| 1 1 I I I 1 1 1 1 1 1 I I I 1 1 1 1 -JT) 

~-&-_ G __e-©--©--G--e-€>' _ 




Al/l 



1 



10 100 1000 
number of initial edges 

scheme I - - e - - scheme II 



1 I 1 — i i i ■ 1 1 ■ | 1 — i i i ■ 1 1 1 1 1 — I I I I llll 

01- G--e-e L -o-^3--e-e--o--G--e 
0.01 - 
0.001 - 
le-04 - 
le-05 - 
le-06 - 
le-07 - 
le-08 - 
le-09 - 

le-10 1 — 




10 



100 



1000 



number of initial edges 
- - + - - scheme III 



Figure 6. Maximal relative errors in Kruskal space-time (time- 
like initial data). Left: values of R~ 2 at the nodes (scheme I) and 
values of a at the space-like edges (schemes II, III). Right: invariant 
lengths of space-like edges. 



6. Discussion 

From the results presented in the last section we conclude that the schemes I and 
III provide very good results. Especially scheme III converges quadratically to the 
analytical solution in every simulation we performed. Also in scheme I most of the 
simulations led to quadratically convergent results. There are only a few regions in 
Kruskal space-time where the errors of this scheme become quite large. This is the 
the case for R w 3M. The reason for this is, that there two solutions of the discrete 
equations are close to each other. The root finding algorithm then sometimes 
chooses the wrong one. However, this should be solvable with an optimised choice 
of starting values for the Newton iteration, but we did not find a practical way to 
obtain such starting values. 

The errors obtained with scheme II are on the other hand quite large. This 
scheme is in some situations not convergent and its behaviour strongly depends 
on the initial data. We conclude that it is not feasible for numerical calculations 
and thus seems to be ill-posed. Now, the question is what the reasons for these 
problems are. 

The difference between the schemes II and III was the direct implementation of 
the Hodge operator in scheme III while scheme II made use of the equations (fT9]l . 
which have been added to the system in order to encode the duality between the 
forms a and (3. It turns out, that far from the horizon in the exterior the values of 
a. on the space- like edges are much smaller than the values of (3, while far from the 
horizon in the interior it is the other way round. Near the horizon the two 1-forms 
interchange their role in this sense, and are hence of comparable size. 

Near the horizon the solution of (fT9|) largely differs from a. = ±/3. This prop- 
erty appeared in all simulation where the initial hypersurface crossed the horizon, 
especially it seems to be independent of the location of the nodes in the initial 
hypersurface. The origin of this behaviour seems to be the ill-conditioning of the 
linear system corresponding to Thus (|19p is a bad implementation of the 

Hodge operator in that region. 

Clearly the Hodge operator (f2Tj) leads to a convergent scheme. However, ([2Tj) 
was inspired by the dual mesh of electrodynamics. In a (1 + l)-dimensional Lorentz 
geometry it is convenient to define the dual of a light-like edge to be the edge itself. 
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Yet, we wanted to avoid the construction of a dual mesh. The first reason for 
this decision is that without a dual mesh the discrete system decouples into small 
systems for every face. The second reason is that the dual of an initial edge does 
not lie in the initial hypersurface and hence we do not know how to specify initial 
values for it. The third reason is that the dual cells are in general not simplices 
anymore, which causes problems in defining a natural discrete exterior product 
without ambiguities. 

The discrete Hodge operator (fT§|) also made it necessary to use the regauging 
described in section 13.41 Although we cannot make definite statmcnts about this 
procedure yet, it seems that it is also a source of errors in scheme II. These prob- 
lems probably arise because the notion of a gauge transformation for the discrete 
equations has not been clarified completely so far. This is planned to be the topic 
of future investiagtions. 

7. Conclusion 

In this article we presented first results of the application of discrete differential 
forms in General Relativity. It was shown, that the method is quite promising. 
Several schemes were found whose results are close to the analytic solution and the 
errors of which converge quadratically to zero. 

We discussed that one has to be careful with the definition of discrete Hodge 
operators and that the notion of discrete gauge transformations is not completely 
understood. Making a wrong decision in these fields can lead to results with big 
errors. 

Now one has to get a better understanding of gauge transformations for the 
discrete equations, a more general definition of the Hodge operator is necessary to 
be able to include matter fields and the method must be applied in space-times 
with smaller symmetry groups, in order to get physically more relevant solutions. 
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Appendix A. Derivation of the reduced system 

A.l. Spherical symmetry. In this appendix we will derive the reduced exterior 
system (0. So we assume the existence of an isometric action of the rotation 
group G = 5*0(3) on the space-time manifold i.e., for each element a € G there is 
an isometry <j> a of (A4 , g) . We assume that the orbits of this action are 2-dimensional 
submanifolds except for the fixed points, which form a 1-dimensional submanifold, 
called the origin O. Clearly, the 2-dimensional orbits are round spheres, they carry 
an induced metric which is a constant multiple of the unit-sphere metric. 

Given a point p E M \ O and the orbit S p through p we can split T V M. into 
the 2-dimensional tangent space H p of S p at p, which we call the horizontal space 
and its orthogonal complement V p , the vertical space@ The isometric action maps 



For a point p G O this decomposition of T P M is not possible, but this is no problem, because 
in the simulations the origin was not included in the computational domain. 
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vertical spaces into vertical ones and horizontal spaces into horizontal ones. Fix 
a point p then the isotropy group of p is isomorphic to SO (2) and the symmetry 
action defines a representation of SO(2) on H p . This is isomorphic to the defining 
representation. On the other hand, the isotropy group of p acts trivially on V p . 
This follows from the fact that such an action defines a homomorphism from 5*0(2) 
to 5*0(1, 1). Due to the different topologies of these groups this must be trivial. 

We will be concerned with invariant objects, i.e., objects which are mapped onto 
themselves by this symmetry action. Consider an invariant function /. It satisfies 
0a/ = f ° 4>a = f for every rotation a e G. This implies that / must be constant 
on each orbit because for any two points on a given orbit there is a rotation which 
maps one to the other. 

An invariant vcctorficld V coincides with each push- forward, i.e., 4> a *V = V. 
Suppose V is horizontal then it follows from the free action of the isotropy groups 
on the horizontal spaces that V must indeed vanish because at each point p it 
must coincide with all its images under elements of the isotropy group. Thus, an 
invariant vectorfield cannot have horizontal components so it is always vertical. It 
follows from this that the covariant derivative of an invariant vectorfield along an 
invariant vectorfield is again invariant, hence vertical. Also the commutator of two 
invariant vectorfields is vertical. This implies that the subbundle of the tangent 
bundle consisting of vertical vectors forms an integrablc distribution: any vertical 
vectorfield U can be written as a linear combination of invariant vectorfields Xi 
and X 2 

(31) U = C^Xi + U 2 X 2 
so that the commutator of two such vectorfields is 

(32) [U, V] = U(F i )X i - V([/ J )X, + U l V k [X,, X fc ] , 

hence it is vertical. The maximal integral manifolds will be called vertical surfaces. 
This discussion shows that the space-time has the topology M — Mi x S 2 where 
Ali is a two-dimensional manifold. We define the two canonical projections 

(33) w.M^Mu p:M^S 2 

mapping on the first, resp. second factor. 

We can set up adapted local coordinates as follows. Fix a point p and assign to 
it the coordinates (x°, x 1 , x 2 , x 3 ) where (x ,^ 1 ) are local coordinates near n(p) e 
Mi and (x 2 ,x 3 ) are local coordinates near p(p) £ S 2 . In these coordinates the 
projections are ir(x° , x 1 , x 2 , x 3 ) = (x°, x 1 ) and p(x°, x 1 , x 2 , x 3 ) — (x 2 ,x 3 ). 

We recall the following facts. The projections ir and p induce isomorphisms 
between V p and T^i p \M\ resp. H p and T p ^S 2 . The group action happens only on 
the second factor, i.e., 7r o </> a = w. This implies that p- forms on M which are pull- 
backs from A^i are invariant and annihilate horizontal vectors. Conversely, a p- form 
on M which is invariant and annihilates horizontal vectors is the pull-back of a form 
on Mi. This is true for any multilinear map. Furthermore, invariant vectorfields 
(which are necessarily vertical) project down to a well-defined vectorfield on Mi 
and each vectorfield on Mi corresponds to an invariant vectorfield on M. 

The metric on M can be written in the form 



(34) 



.9 = .91 + .92 
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where gi (#2) is non-zero only for vertical (horizontal) vectors. This shows that the 
topological decomposition is also an orthogonal decomposition. Furthermore, for 
any infinitesimal isometry £ we have 

(35) = C$g = C^gi + C^g-2 

and the orthogonality properties imply that the two terms on the right have to 
vanish separately, 

(36) C^gi = = C i g 2 . 

Now, the facts that C^gi — and that gi(X, •) vanishes for any horizontal vector 
X imply that g\ is the pull-back of a metric h on M\. 

The metric 172 is conformal to the metric S on the unit sphere S" 2 where the 
conformal factor may depend on the coordinates (x®,x l ). Thus, we can write the 
space-time metric in the form of a warped product 

(37) g = Tr*h + R 2 p*S where R 2 : Mi — ► K 

We can now set up an adapted tetrad (0 l )i=o:3- To this end we choose a frame 
(tf ,!? 1 ) on (Mi,h) and a frame (tf 2 ,# 3 ) on {S 2 ,S) and set 

(38) (0°, 1 , 6 2 , 8 3 ) = (vr*i9°, Tr*tf\ Rp*tf 2 ,Rp*'d 3 ). 

Using this tetrad in the first structure equation one finds after some calculation 
that the connection forms are 

w°i=7r*ta>, u; 2 o -/ o 2 , u: 3 = foO 3 , 

u 2 3 =p*zu, w 2 1 =/ 1 2 , u> 3 i = f 1 9 3 . 

Here, uj resp. vj are the connection forms of the metrics on .Mi resp. S 2 and 

(40) a = f 6° + fxd 1 = dR/R 

is an invariant 1-form on Mi. 

Appendix B. Reduction to 1+1 dimensions 

Our goal is to express the field equations (JSJ) as equations on .Mi. The easiest 
way to achieve this goal is to compute the Nester-Witten and Sparling forms from 
the connection forms (|39|) . This yields the following result 

(41a) L = 2f 1 6 2 9 3 -u: 2 3 e 1 , L 2 = -fi0°0 3 - f l 6 3 - 

(41b) L x = 2f e 2 9 3 + u; 2 3 0°, L 3 = /i0°0 2 + / o 1 2 +cA0 2 

for the Nester-Witten form, while the Sparling 3-form is given by 

(42a) S = 2/ o /i0°0 2 3 + (f 2 + f 2 )0 1 2 O 3 + 2f O Lj° 1 2 3 - u a 3 w O i0 , 

(42b) Si = 2/ o /i0 1 2 3 + (/ 2 + / 2 )0°0 2 3 + 2f 1 u>° 1 2 3 + wW 1 , 

(42c) s 2 = hu\o°e 3 - fo^^e 2 + h^ie'e 3 - fr^^e 2 + u>\u 2 3 e 2 , 

(42d) S 3 = -/ o w o i0°0 2 - f u) 2 3 1 8 3 - fi^^O 2 - /ia; 2 30 o 3 + lj Q iu 2 3 3 . 

Let us now compute the exterior derivative of Lq. We obtain 

(43) dL = 2d/i0 2 3 + 2/id(0 2 3 ) - dw 2 3 1 + w 2 3 d0 1 . 
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Since 2 8 3 = R 2 p*(tf 2 -d 3 ) we find 

(44) d(0 2 3 ) = ^8 2 3 

R 2 

and for dw 2 3 we obtain 

(45) dw 2 3 = p*dvj = p*fl 

where f2 is the curvature 2-form of the unit sphere which is given by n = iTtf 3 . 
Hence, 

(46) dw 2 3 = p*(tf 2 tf 3 ) = -^0 2 O 3 . 

R 2 

Furthermore, making use of the structure of the connection forms we have 

(47) dO 1 = -uj^e - u) 1 2 2 - uj 1 ^ 3 = -^00°. 

This shows that with 1 also dO 1 is invariant (this also follows from the fact that 
6 1 = tt*^ 1 so that dO 1 = Ti^dtf 1 ). Taken together, we have 



(48) 



dL = (2dh + - -^e 1 ^ e 2 e 3 + u 2 ^ 1 

= (2foh0° + (f 2 + fDO 1 + 2/ouA) e 2 e 3 - u, 



Using (|4"5)) and (|42ap we write the equation dL^ = So (which is (|5b[) for i = 0) in 
the form 

2d/ a + 2.a^ - ^ 2/oAe - (f 2 + fDe 1 - 2f u\\ e 2 e 3 = o 



R 2 R 2 

and contract with two horizontal vectors. Then we find the equation 
(49) 2d/ x + 2/^ - -^O 1 - 2foh9° - (J 2 + fDO 1 - 2/ocA = 0. 



Conversely, if (|49|) holds then the field equation d_Lo = So is satisfied. In a similar 
way we can treat the equation dii = Si, i.e. (|5b|) for i = 1. We obtain 



(50) 2d/ + 2/ ^ + -^9° - 2/ohe 1 - {f 2 + fl)0° - 2fru>\ = 0. 

For i — 2,3 equation (|5bj) has to be treated differently. Let us consider the case 
i = 2. We first define the invariant 1-form 

(51) p = he 1 + he 

and then we derive from (|41a|) 

(52) dL 2 = -dip + w°i)0 3 + ((3 + cA)d0 3 = -auAfl 3 -{13 + u)°i)u) 3 2 2 . 
Using the fact that 

(53) d0 3 = -Lj 3 Q e° - a; 3 !© 1 - lj 3 2 2 = -u 3 2 6 2 + a.6 3 
we can write the equation dL 2 = S 2 as 

(d/3 + dcA +ol0) 6 3 = 0. 
Contraction with one horizontal vector shows that this equation implies 

(54) d/3 + du> O i+a/3 = 
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and, conversely, if (|54|) is satisfied then also dL2 = S2 holds. The equation dL% = 
S3 does not contribute anything new. 

Now we can collect the equations. The contents of the first structure equation 
is the specific form of the connection forms, the first structure equation on S 2 and 
the first structure equation on Ai\. Since all the relevant quantities are invariant 
they are pull-backs from .Mi. So the equations really live on Mi. However, in 
order not to complicate the notation we will continue to write them as equations 
on M. , keeping in mind that everything has to be regarded as a pull-back from the 
2-dimensional manifold M.\. Then we have 

de° + u° 1 e 1 = 

( 55 ) 1 1 n 

Furthermore, we have the relationship between the 1-form a. and the differential 
of R 2 

(56) dR 2 = 2aR 2 . 

That means the integral of a. along a curve is the same as the difference of the 
values of (logi?) at the boundary of that curve. Hence a is related the velocity 
with that the area of the spheres changes. 

Then the field equations written in terms of fa, fi and /3 = JqO 1 + f\0° are 

(57) d/ - w/i + ha + 1 (y 2 - f 2 + ^\e° = 0, 

(58) d/i - «/„ +h<*~\ (fo - fi + ) Q 1 = 0, 

(59) d/3 + du; + a/3 = 0. 

Note, that contracting the first of these equations with 6° and the second with 
9 1 yields the integrability condition da = 0. On the other hand, using these two 
equations to compute the differential of f3 we find 

(60) d/3 = -2af3 - -^0 V 

R z 

and, therefore, 

(61) duj = af3+^-e e 1 . 

Rr 

Note, that a and j3 are not independent. They contain the same information. In 
fact, we have j3 = -ka. This relationship can be expressed in the present case also 
as two 2-form equations 

(62) aO 1 + [30° = 0, a0° + [30 1 = 0. 
This concludes the derivation of the reduced equations. 
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